Vortex Pinning and the Non-Hermitian Mott Transition 
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The boson Hubbard model has been extensively studied as a model of the zero temperature 
superfluid/insulator transition in ''He on periodic substrates. It can also serve as a model for vortex 
lines in superconductors with a magnetic field parallel to a periodic array of columnar pins, due to 
a formal analogy between the vortex lines and the statistical mechanics of quantum bosons. When 
the magnetic field has a component perpendicular to the pins, this analogy yields a non-Hermitian 
boson Hubbard model. At integer filling, we find that for small transverse fields, the insulating 
phase is preserved, and the transverse field is exponentially screened away from the boundaries of the 
superconductor. At larger transverse fields, a "superfiuid" phase of tilted, entangled vortices appears. 
The universality class of the transition is found to be that of vortex lines entering the Meissner phase 
at Hci, with the additional feature that the direction of the tilted vortices at the transition bears 
a nontrivial relationship to the direction of the applied magnetic field. The properties of the Mott 
insulator and fiux liquid phases with tilt are also discussed. 



I. INTRODUCTION 



The importance of interactions for vortex lines in su- 
perconductors has been well known since Abrikosov's the- 
ory of the mixed phase of type II superconductors. The 
vortices form a regular triangular lattice due to their re- 
pulsive interactions. With the discovery of the cuprate 
high-temperature superconductors, however, it became 
clear that thermal fluctuations can play an important 
role in the behavior of the flux lines arrays. ||l| The com- 
petition between interactions, pinning, and thermal fluc- 
tuations gives rise to a wide range of novel phenomena. 
For example, it is now clear that the triangular flux lat- 
tice melts via a first order transition Plin clean materials, 
giving rise to an entangled flux liquid |3|] over a significant 
part of the phase diagram. 

Here we investigate the behavior of an interacting 
system of flux lines in a periodic array of columnar 
defects. Traditionally, columnar defects are made in 
high-temperature superconductors by heavy-ion bom- 
bardment. The ions pass through the sample, leav- 
ing amorphous damage tracks in their wake, which act 
very effectively as pinning sites for the flux lines, provided 
the lines are aligned with the defects. Pinning is efficient 
in this case because the dimensions of the flux lines are 
comparable to those of the pins: both the vortices and the 
columnar defects are extended over the length of the sam- 
ple in one direction, and the vortex core size (~ 20 A) is 
comparable to the width of the damage tracks (~ 60 A). 
By covering the sample with a mask with a periodic ar- 
ray of holes through which the ions can pass, a periodic 
array of pins (or clusters of pins) could be formed. 

There are, however, other promising approaches to 
forming a periodic array of columnar defects in which the 
artificial pinning centers are able to accommodate mul- 
tiple vortex lines per pin. Such structures have in fact 
already been made both in thick |^ and thin samples, ||l^] 
and one goal of this paper is to understand effects of ther- 
mal fluctuations and entanglement when the materials 



are very thick. The work of Cooley et al. ^ is of partic- 
ular interest to us, because their samples can potentially 
be used in two different ways to test our theories. Their 
samples have a regular lattice of pins within an "island." 
The islands are themselves spaced periodically through 
the sample. We may use our theory to predict behav- 
ior of the pins within a single island; or, we may view 
each island as a large pinning center capable of accom- 
modating many vortices and apply our theories to that 
as well. The lattices are triangular in either case. We 
here describe effects expected to be observed in square 
lattices of pinning sites; it is straightforward but tedious 
to generalize our results to triangular lattices as well. 

Another motivation for examining flux lines in the 
presence of an ordered array of columnar pins is to bet- 
ter understand computer simulations such as those by Li 
and Teitel jl^] and others which attempt to model 
flux lines in cuprate superconductors. These authors use 
an underlying discretization (such as a stack of triangu- 
lar lattices) for the simulations. Although layering along 
the c-axis is an important feature of real cuprates, the 
relatively coarse mesh of allowed sites for vortices in each 
layer is clearly an artifact of the simulations. As we shall 
see, the presence of such periodicity can cause significant 
changes (such as a transverse Meissner effect) in the be- 
havior of the flux lines. The very existence of a periodic 
"Mott insulator" phase in these simulations (see below) 
with its massive phonon modes and infinite tilt modulus, 
is an artifact of the lattice discretization. In the absence 
of a lattice, the Mott insulator phase would be replaced 
by an Abrikosov crystalline phase with finite elastic con- 
stants. Our investigation can begin to delimit the circum- 
stances under which lattice approximations fail as models 
of real materials, and to predict what kind of behavior 
might be observed when they do fail. 

As a third motivation, we hope that studies such as 
this will lead to a better understanding of interacting 
flux lines with tilted field in a random array of colum- 
nar defects. IfqjlJ] Although the extension of the stan- 
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dard theory of disordered superfluids to vortex lines 
with tilted magnetic fields is straightforward in princi- 
ple there are still poorly understood aspects of 
this problem. 

To understand how our model differs from conven- 
tional lattice boson Hubbard models, |jl^ consider the 
transfer matrix generated by the classical partition func- 
tion for wiggling vortex lines in d -I- 1 dimensions. This 
transfer matrix is the exponential of a d-dimensional non- 
Hermitian quantum Hubbard Hamiltonian 

The connection between this model and the physics of 
flux lines will be reviewed in more detail in Sec. Q. Here, 
we simply note that the {a\,aj} are boson (i.e., vortex) 
creation and destruction operators, the are unit vec- 
tors connecting nearest neighbor lattice sites on a square 
lattice representing the positions of columnar pins in d 
dimensions, and the first sum is over neighboring pairs of 
pins. The parameter t controls vortex hopping between 
columns, /i — /i(-ffz) is a chemical potential that depends 
on the component of the external field Hz parallel to the 
columns, and U represents the penalty for multiple vor- 
tices on a single defect. The external tipping field Hj^ 
determines the non-Hermitian hopping asymmetry h via 
h = (/!)oHj^ao/(47rA:BT), where ao is the lattice constant. 
In the presence of this non-Hermitian hopping, which 
would be impossible in conventional many body quantum 
mechanics, the excitation spectrum can become complex. 
The complex spectra which arise when J7 = and /i varies 
randomly from site to site has aroused considerable in- 
terest recently. In this paper, we study how this 
non-Hermiticity affects excitations and other features of 
the Mott insulator and superfluid phases. 

Before discussing the interacting case, consider a sys- 
tem of noninteracting flux lines in the presence of a pe- 
riodic array of columnar pins. At very low temperatures 
in thin samples, each vortex will be tightly bound to a 
particular pin, and fluctuations away from the pinning 
sites will be small. However, a flux line can gain entropy 
by hopping from one defect to another. For an infinitely 
thick sample, the vortices will in fact be delocalized at 
any nonzero temperature. In the quantum analogy, this 
delocalization is a consequence of the periodicity of the 
pinning potential and Bloch's theorem. Although the 
distance in the z-direction between intercolumn hops di- 
verges as the temperature is lowered to zero, p9[| we will 
assume sufficiently thick samples and high temperatures 
so that an individual vortex hops repeatedly as it tra- 
verses the sample. Such delocalized states are easy to 
tilt and will have a finite linear resistance, because trans- 
port in the presence of a Lorentz force will be aided by 
finite concentration of kinks in equilibrium. p3| Hopping 



from one columnar pin to another cannot be neglected at 
finite temperatures in thick samples, even if such events 
are relatively rare. 

We now discuss qualitatively the interacting case. Con- 
sider an external magnetic field Hz such that there is ex- 
actly one flux line per columnar pin, i.e., Bz = 4'o/o,q- 
If the pins are strong enough, the ground state will be a 
square lattice of vortices rather than the usual triangu- 
lar array that is observed in the absence of pins. Ther- 
mal fluctuations allow the lines to make short excursions 
away from the columnar pins which lowers the effective 
strength of the pins, jl^ More substantive changes arise 
if thermal disorder also leads to hops from one pin to 
another, as for the noninteracting flux lines discussed 
above. However, the interactions inhibit hopping, due 
to the extra energy cost of double occupancy. Q Al- 
though hopping increases the entropy per unit length, 
this is insufficient to overcome the energy barrier at low 
temperatures. Above a certain critical temperature Tc, 
however, the entropy gain dominates, and the flux lines 
can hop freely from pin to pin. In the quantum analogy, 
this phase change is the boson Mott transition be- 
tween a 2c? Mott insulator (square vortex lattice) and a 
2d superfluid (entangled vortex liquid) . The temperature 
T for the vortices plays the role of h for the bosons. (See 
Sec. ig.) 




FIG. 1. An example of a correlated hopping event. Such 
hops delocalize flux lines, but cannot mediate vortex 
transport or tilt. 

Correlated hopping events such as are shown in Fig. 1 
can lower the entropy without forcing two flux lines to 
occupy the same column, and always delocalize lines in 
sufficiently thick samples even in the crystalline phase. 
However, such defects cannot mediate vortex transport 
or tilt — if one flux line moves to the right, another must 
move to the left, so there is no net change. Thus it will 
not be until the temperature rises above the critical tem- 
perature mentioned above that there will be any linear 
resistance or response to a perpendicular magnetic field. 

Consider now raising the applied magnetic field H = 
HzZ, always keeping it parallel to the column direction. 
In the absence of columnar pins, the lattice spacing would 
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simply shrink; however, if the pins are strong enough, a 
change in the lattice constant would sacrifice the energy 
benefit of commensurability with the lattice of pinning 
sites. An alternative response is the introduction of ex- 
tra vortices into the commensurate lattice. Depending 
on the pinning strength, the additional vortices will ei- 
ther doubly occupy certain columnar pins, or else occupy 
interstitial sites; either possibility costs an extra energy 
cost per unit length. We assume for simplicity that the 
cost of double occupancy is lower |2^, and call the re- 
sulting extra vortex a "particle" excitation. This "par- 
ticle" has an energy per unit length of order the cost 
of double occupancy; however, it can gain entropy by 
hopping from pin to pin, and thus lower its free energy. 
In fact, its free energy can become negative. Follow- 
ing the treatment of interstitials and vacancies in a con- 
ventional Abrikosov lattice in Ref. an estimate of 
the boundary above which particles are favored arises 
from approximating the partition function of a single ex- 
tra vortex by 4^/'^ exp(— epL/T), where ep — ep{Hz) is 
the energy cost per unit length of a particle, L is the 
length in the z direction, and 1^ is the distance between 
hops of the particle from one site to another. If Ep is 
the extra energy of an particle in zero external magnetic 
field H, then ep{H) = Ep - <j)oH^/47r. The factor of 4 
comes from the four directions available for hopping on 
a square lattice. Thus the free energy of the particle will 
he F fii €pL — T{L /Iz) In 4, which becomes negative above 
a temperature 



In 4 



(1.2) 



Above this temperature, particles are favorable; below 
it, however, the system cannot respond to a change in 
and the extra magnetic field is screened out. Sim- 
ilar screening arises if one attempts to produce vacan- 
cies or "holes" by reducing the magnetic field, with 
Fh « ehL—T{L/lz) In 4, where eh is the hole energy in the 
absence of an external magnetic field. Thus for a range of 
applied fields Hz , the actual field (as measured by the 
vortex density) will be locked in at the value Bz — (po/a^, 
much as in the Meissner effect where the field is locked 
in at _B = for a range of applied Hz ■ Similar arguments 
show that if we apply a small field Hj^ transverse to the 
columns, isolated flux lines will not be able to tilt in re- 
sponse to it, and the field will again be screened out, i.e., 
Bi = 0. II 

In the rest of the paper, we examine the phase tran- 
sitions described qualitatively above. The phase dia- 
gram is shown in Fig. 2. The "commensurate" Meissner- 
like phase discussed above for vortex matter is called 
a "Mott insulator" in the context of quantum bosons. 
With Hj^ = 0, this problem has been studied by Fisher 
et al. in the context of real quantum bosons on a lat- 
tice. In the presence of a transverse field, however, the 
formal analogy between the statistical mechanics of vor- 
tex matter and (2-|-l)-dimensional quantum mechanics 
leads to the non-Hermitian model of Eq. (^) . 
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FIG. 2. (a) Schematic phase diagram for = 0. The 
circles at the tips of the Mott lobes are the multicritical 
points, (b) Schematic phase diagram at fixed temperature 
T below the multicritical points in (a). 

In Sec. H we describe in more detail our model of in- 
teracting vortex lines. We study the Mott insulator to 
fiux liquid phase transition induced by a perpendicular 
field in Sec. [II, using mean field theory and a more so- 
phisticated renormalization group treatment. We argue 
that in the presence of a transverse field, the transition 
is characterized by the penetration of "tilted defects," 
which are particles or holes with net average motion in 
response to the transverse field, as shown in Figs. 3 and 



4. More support is given for this picture in Sec. IV, 
where the low temperature Mott insulator phase is ex- 
amined. We determine excitation gaps for hopping in 




FIG. 
feet. 



3. An example of a right-moving tilted particle de- 
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the presence of the non-Hermitian perturbation. We also 
calculate a quantity analogous to the London penetra- 
tion depth in the Meissner phase, which is the length 
over which screening of the transverse magnetic field be- 
comes effective near the boundaries over the sample. In 
Sec. 0, we examine properties of the superfluid phase, 
particularly the way we expect the magnetic field to pen- 
etrate the sample near the transition, and the excitation 
spectrum. 




FIG. 4. An example of a right-moving tilted hole defect. 



II. MODEL 

We begin with a model free energy for N vortex lines 
in a sample of thickness L in the z direction (perpendic- 
ular to the Cu02 planes) in the presence of a periodic 
square array of columnar defects aligned in the z direc- 
tion with nearest neighbor separation a. We allow an ex- 
ternal magnetic field Hj^ transverse to the defects as well 
as the usual parallel field H|| . The energy then reads 



dz 



dz 



IJ2 fv{\v,{z)^v,{z)\)dz 
/ VD[v^.{z)]dz 



An 



with 



M 



yz5[r^(z)]=^Fi[|r^(z)-R.|]. 



(2.1) 



(2.2) 



Greek letters denote the flux lines while Roman letters 
denote the columnar pins. V{r) is the (repulsive) inter- 
action potential between flux lines, which can be taken 
local in z, while Vi describes the (attractive) interaction 
between a flux line and a columnar defect. The param- 
eter ei in the kinetic energylike term arises from the tilt 
energy of the lines. This term is a quantitatively correct 
approximation to the tilt energy provided 



rfr^(-z) 



dz 



(2.3) 



where Mz/M±^ ^ 1 is the mass anisotropy in the un- 
derlying Ginzburg-Landau theory. That same criterion 
insures that an interaction local in 2; is a good approxi- 
mation, jl^ 

There is a useful formal correspondence between this 
flux line problem and flctitious quantum-mechanical 
bosons in two dimensions. [^JT^. The temperature T 
plays the role of Planck's constant ?i, the bending en- 
ergy ei plays the role of the boson mass m, and the 
length of the sample L corresponds to jSTi for the bosons. 
(See Table I for a summary.) However, while real bosons 
have periodic boundary conditions in imaginary time, free 
boundary conditions are more appropriate to the flux line 
system. As we take the thermodynamic limit L 00 
(corresponding to the zero temperature limit for the fic- 
titious bosons), the partition function is dominated by 
the ground state, and the statistical mechanics is insen- 
sitive to the boundary conditions. Even for finite L, the 
only states which contribute to the partition function of 
the vortex lines are symmetrical boson states . W e de- 
fer consideration of finite size effects to Sec. IV B, while 



in the rest of this paper we use periodic boundary condi- 
tions for convenience and take the limit L — s- cx). 



Vortex lines 


Bosons 


ei 


m 




h 


L, 




i?^0o/47r - ei 


M 


Hj_0oao/47rfcsr 


h 


Meissner-likc commensurate phase 


Mott insulator 


Flux liquid 


Superfiuid 



TABLE I. Detailed correspondence of the parameters of 
the flux line system with the parameters of the two- 
dimensional boson system. 

As emphasized in Ref. [Q, the transverse field Hj^ 
plays the role of an imaginary vector potential in the 
boson system. Thus, understanding the fiux line sys- 
tem in the presence of a transverse field is equivalent to 
understanding the quantum mechanics of non-Hermitian 
two-dimensional bosons with interactions and a periodic 
potential. 
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In the limit of a very stro ng columnar pinning poten- 
tial, we can replace Eq. (2.1) by a simplified tight binding 
model for the vortex lines [ [I3|jl^j20| ] such that, written in 
terms of the fictitious boson operators, the grand canon- 
ical partition function is 



Z = Tr 



(2.4) 



where Ji given by Eq. (O) and L is the sample thick- 
ness. The non-Hermitian terms in the effective Hamil- 
tonian reflect a preference for the bosons to hop in the 
direction of the transverse field. We have neglected inter- 
actions between flux lines that are not on the same site. 
This should give reasonable results near filling fraction 1; 
however, if we wanted to describe filling fraction 1/2, for 
example, where we might expect the ground state to be 
a checkerboard arrangement of the flux lines, we would 
need to include further neighbor interactions. Further 
neighbor interactions would also be required to treat the 
lattice models of melting studied numerically in Refs. |l|] 
and Q. 

The parameters of this model, expressed in the original 
flux line language are |T^ , [l^ 







(2.5) 


h = 


Hj_0oao 
47r/eBT ' 


(2.6) 


[/« 


2eoln(F,2/i?), 


(2.7) 


t K. 




(2.8) 



with Vo the strength of the pin and Eij — \/2eiVbao the 
WKB exponent for tunneling from one pin to its neigh- 
bor. The Hubbard U term in Eq. (1.1) contributes an 
energy ^Un{n— 1) when n flux lines are on a single site. 
It represents the difference between the energy of a sin- 
gle n-fold vortex quantum (~ n'^) and the energy of n 
well-separated unit vortices. In this tight-binding limit, 
a Mott insulator phase exists not only at filling fraction 1 , 
but at any integer filling fraction. We will consider Mott 
insulator phases with general filling fraction n. Integer 
filling fraction with n > 2 would arise as an approxima- 
tion to samples irradiated through a mask of holes large 
enough to allow multiple columnar defects at every lat- 
tice site. 



parameter t. As discussed in Sec. |[ the "Mott insulator" 
corresponds to a Meissner-like vortex phase, which resists 
changes in the number of vortices per site. As is evident 
from Eq. ( ^.8[) , small t corresponds to low temperatures 
in the vortex system. Large hopping leads eventually to 
a 2nd order phase transition to a superfluid phase. This 
"superfluid" is a lattice analogue of an entangled flux liq- 
uid for vortex matter. In this section, we first mention 
a simple mean field theory due to Sheshadri, et al. |2|] 
The detailed application of this mean field theory to the 
non-Hermitian boson Hubbard model is described in Ap- 
pendix |A[ We then use a more systematic approach to 
study this transition which will allow us to determine the 
universality class of the transition. We defer a detailed 
discussion of the phases themselves until later sections. 

A. Mean field theory 

One simple way of accessing the mean field theory as- 
sociated with Eq. (^^) is with an ansatz that decouples 
the hopping term to give a single-site Hamiltonian via 
the replacement [E2[ 



(3.1) 



In the remainder of this section, however, we follow a 
different approach that allows a more convenient treat- 
ment of fluctuation effects. The decoupling method, how- 
ever, is conceptually and computationally simpler, so we 
present it in Appendix 



B. Hubbard-Stratanovich transformation 

We here transform the tight binding partition function 
of Eq. (2.4) and Eq. (1.1) into a field theoretic form more 



convenient for calculation following a method used by 
Fisher, et al. |jl^, and earlier by Doniach |^^, to treat a 
quantum model of coupled arrays of Josephson junctions. 
We begin by separating the Hamiltonian Hjy of Eq. (1.1) 
into the single site piece Tio and the off-site hopping term 
Hi: n ^Ho + Hi, with 



U 



hi{ni - 1) - nfii 



i 

Hi = -- '^[Jijaldj + Jjio'^CLi], 



(3.2) 
(3.3) 



III. TRANSITION AND RENORMALIZATION 
GROUP 

Unless fi/U is exactly an integer [corresponding to the 
points where two Mott lobes are tangent in Fig. 2(a)], it 
is expected JlGf that the equivalent quantum system will 
be in a Mott insulator phase with a well defined integer 
site occupation number n for sufficiently weak hopping 



where fii = ajoi is the number operator and Jy is the 
hopping matrix element between sites i and j: 



Jr. 



^gh eij j£ j nearest neighbors 
otherwise. 



(3.4) 



The generalization of the transformations of Refs. |2^] 
and [ p^ to the non-Hermitian case presents no difficul- 
ties. Their effect is to allow the grand canonical partition 
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function of Eq. (2.4) to be rewritten as a functional inte- 
gral over a set of complex functions {ipiir)}, 

2 = l[vMT)Vr^{r)exp[-Sii;)], (3.5) 



with an effective action 



^Tt- exp 









+ ^*(T)a.(r)} 



(3.6) 



J / 



where Tr is the imaginary-time ordering operator, 

Zo Tr|e"''^^"| , (3.7) 



(3.8) 



and the imaginary time dependence of operators is given 
by, e.g., 



(3.9) 



We can see the equivalence of Eq. (3.5) with Eq. (2.4) by 



integrating out the V' fields to obtain 



Z = Zq l Tr exp 



'(3 / dTUiir) 
'o 



(3.10) 



J / 



whic h is just the familiar interaction representation of 
Eq. (O). 



C. Universality classes 

Since (ipi) is an order parameter for superfluid- 
ity p3| , p^ we can study the tran sition to a state with 
(ipi) ^ by expanding Eq. (3.6) in powers of ip- For 
now, we keep only the most relevant terms in a long- 
wavelength approximation, and obtain 

/(P \iduj 
-j^-^\ij(k,uj)\^[r - iCf,uj + ich • k 

+ c'uj^ + D^jhkj]+u J (frdT\^P{r,T)\^, (3.11) 

as shown in Appendix This formula is applicable 
to a Mott insulator phase with n particles per site (i.e. 
n — 1 < fi/U < n). The values of the parameters r(/x, T), 
c', M, and the matrix Dij are given in Appendix ^ while 



p 



El 



Ch 



where 



ao X sinh hx + y sinh hy 
2t [cosh /la; + cosh/iy]2 



Ep = Un — /i. 



(3.12) 
(3.13) 



(3.14) 

Eh= ^i-U{n~l) (3.15) 

are (respectively) the energies of "particles" and "holes" 
superimposed on a Mott insulator state with a fixed value 
of n, in the absence of hopping. 

Note that and Ch can be positive, negative, or zero. 
In the case examined in Fisher et al, there is no trans- 
verse field, i.e., Ch = 0. If = as well, then the uni- 
versality class of the transition which occurs for r « is 
that of the XY model in three dimensions, with nontriv- 
ial critical exponents due to fluctuations. When > 0, 
the transition at r w is characterized instead by mean 
field exponents with logarithmic corrections. |1^,|2^ This 
is also the universality class which describes the penetra- 
tion of thermally excited vortex lines from the Meissner 
phase near Hd ■ ^ Near this transition, extra vortex lines 
enter into the sample and behave (at long length scales) 
like a dilute gas of weakly interacting bosons. When 
< 0, the universality class is again like flux penetra- 
tion near Hd, with the magnetic field reversed. Now, 
holes enter the sample and behave (at long wavelengths) 
like a dilute system of linelike excitations in the negative 
z direction. 

When a tran sverse field is applied, Ch ^ 0, and we see 
from Eq. (3.11) that the transition again resembles the 
physics of dilute linelike excitations near Hd, but in a 
direction given by the three dimensional vector 

V = Ch + c^z. (3.16) 

When r(/i, T) < 0, tilted holes or extra vortices enter 



the sample, at an angle determined by Eq. (3.16). See 
Figs. 3 and 4 for illustrations of tilted particle and hole 
defects, and see Fig. 5 for the range of parameter space 
in which we expect each type of defect. The renormal- 
ization group calculation that describes the transition, 
including effects of the underlying square columnar lat- 
tice of preferred sites, is given in Appendix 

We confirm this interpretation of the physics when 
we examine the Mott insulator and superfluid phases 
in Sees. |^ and 0. In the Mott insulator, we calcu- 
late the energies per unit length of defects constrained 
to have a specific average direction, and find that only 
those oriented in direction v have vanishing energies per 
unit length as we approach the transition. Thus, near 
the transition, only defects in this specific orientation 
will have favorable energy and proliferate. In the super- 
fluid phase, we calculate the magnetic field change 5'B 
in excess of the field Bq locked into the Mott insulator 
phase, and find that it lies in the v direction as well, cor- 
roborating our assertion that the defects penetrate at an 
average angle 
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9,, = tan 



tan 



1 / |Ch| 



(3.17) 
(3.18) 



where Oy and (pv are polar angles with respect to the z 
axis. 



left-moving 
particles 



right-moving 
holes 



right-moving 
particles 



left-moving 
holes 



FIG. 5. The range of parameter space in which we expect 
to see various kinds of defects as we cross to the superfluid 
side of the Mott insulator/superfluid transition. T he pa - 
rameters and Ch are defined in Eqs. (|3.12|) and (3.13). 



tan 



kBT\h\_ 



(3.22) 



for the polar angles, where 



Afj.0 = - [Eh - Ep + 2i(cosh/i:r + coshhy)] (3.23) 

(3.24) 



M - Mo, 



with 



^0 = U{n — 1/2) — t{coshhx 



cosh hy). 



(3.25) 



From the equations for the azimuthal angles, we can see 
that the defects will tend to penetrate along the crystal- 
lographic axes defined by the square lattice of columnar 
defects. This is because the kink energy is smaller for 
hopping along the nearest neighbor directions, and is fur- 
ther emphasized by Fig. 6, which shows the relationship 
between (f>y and 0h- To better compare 9y and On, we 
note that at the transition, we can write the ratio 



tan 9y 
tan 9 



H 



altU 



/(A^,h), 



(3.26) 



where 



sinh hx + sinh hy 



(3.27) 



D. Relationship between defect angle and the 
direction of the applied magnetic field 

Note that these angles are not the same as the di- 
rection of the excess applied field, defined to be the 
additional field beyond that which produces the Mott 
insulator background. We denote this excess field by 
SH = H — i/oz. Because there is actually a range of 
fields that produce the Mott insulator phase (see Fig. 2), 
we define Hq so that = when SHz = 0. Then a com- 
parison of the direction of the defects with the direction 
of SH near the transition r = gives 



tanffl,, 



sinh h 



V 



sinh hr^ 
h„ 



tan 0H = 7^ 



for the azimuthal angles, and 
aot{U + ijl) 



(3.19) 
(3.20) 



tan^,, — 



sinh /ij; + sinh hy, (3.21) 



and fi and h are related by the constraint that we are at 
the transition r = 0. The prefactor a^tU / [kBTY gives 
the order of magnitude of the ratio of tan 9y to tan 6'h, as 
the (dimensionless) function /(/i, h) is of order 1 unless 
Ihl » 1. 




45 90 135 180 225 270 315 360 

FIG. 6. A plot of the azimuthal angle of the average di- 
rection of defects versus the azimuthal angle of the excess 
applied field, for a fixed |Hx|. 
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E. Scaling analysis near the multicritical point 

A formal scaling analysis determines the shape of the 
transition line near the multicritical XY transition. Gen- 
eralized to Ch 7^ 0, the scaling form for the singular 
part of the free energy fs near the multicritical point 
c = Ch = 0,(5 = r — rc = reads 

fs{6,Ch,c) ^ 6-'^^(^6l/''--,Ch6^^c^6^-) (3.28) 

for a length rescaling parameter b. Ah and are the 
renormalization group eigenvalues for Ch and respec- 
tively. Although the scaling prefactor is in 
general, we use the fact that z = 1 near the XY transi- 
tion. 

If we choose b such that — 1, we obtain 



(3.29) 



where $±(x, z) is a scaling function that depends on 
whether (5 > or (5 < 0. 

A straightforward generalization of arguments in 
Ref. yields Ah = A^ = 1 and Eq. ( 3.29| ) becomes 



/.(<5,Ch,c^)~ |<5|'^''--<i>±(ch|<5| 
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■). (3.30) 



The function <&-(x, z) is singular along a curve [xc(2;), z\, 
which is necessary to describe the superfluid/Mott insu- 
lator transition away from the special multicritical pint 



at Cn 



Ch = 0. Thus Ch 



on the transi- 



tion line. Since vxy < 1 in three dimensions, there is no 
cusp in the phase diagram, in contrast to the Bose Glass 
phase. 



IV. MOTT INSULATOR PHASE 



J-i(k) 



1 



1 



2t cos(fca;ao -I- ihx) + cos{kyao + ihy) 



(4.2) 



for a square lattice of colu mnar pinnin g sit es, and Ep 
and Eh are defined by Eqs. ( 3.14 ) and ( 3.15 ). T he p ar- 
ticle energy Ep is the energy arising from Eq. (1.1) of 
a "frozen" extra vortex, i.e., with no contribution from 
the entropy gain due to wandering of the particle. In 
the quantum analogy, this entropy is represented by zero 
point motion. Similarly, Eh is the hole energy associated 
with a missing vortex frozen at a single site of a Mott 
insulator state with n vortices per site. 

We remark that the quadratic expansion in ijj is suf- 
ficient to determine the Mott insulator/superfluid phase 
boundary within mean field theory. The Mott insulator 
phase is stable provided the coefficient of the |i/;(k, lj)P 
term is positive for all k and oj. Thus the mean field 
transition is determined by the first mode that vanishes, 
which is the k = 0,w = mode. Therefore, the transi- 
tion occurs a t r(/i , Hj^, T) = 0, where r, the parameter 
given in Eq. ( [3.11 ), is given by 



1 



1 



2t cosh h-r + cosh h„ 



n 



1 



E,, 



(4.3) 



The decoupling method also leads to a transition at 
r = 0, as shown in Appendix 



A. Energy gaps for particles and holes 

Because quantities deep in the Mott insulator phase are 
well described by Gaussian fluctuations, we focus on the 
quadratic approximation, and describe later how the re- 
sults change close to the Mott insulator/superfluid tran- 
sition. 



The renormalization group analysis described in Ap- 
pendix ^ allows us to relate quantities near the transition 
to quantities deep within the Mott insulator or superfluid 
phases. The Mott insulator is characterized by the en- 
ergy gaps for excitations. The magnitude of these gaps 
will control tra nsport properties of the sample, as exem- 
plifled in Sec. IV A3 where we calculate the nonlinear 
resistivity in the Mott insulator phase. 

We can study these gaps far from the transition by ex- 
panding Eq. (3.6) to lowest order in ip, keeping all terms 
quadratic in ^, not just the most relevant in a small k 
and uj expansion. As shown in Appendix ^ 



5(V')« ^o(^) = Pal 
X J^i(k)- 



d^kdcj 

n 



lV'(k,c.)p 

n -I- 1 



Eh + iuikBT Ep — iuikBT 



(4.1) 



Here (J~i)(k) is the Fourier transform of (J given 
by 



1. Gaussian fluctuation effects 

To determine the free energy cost per unit length of 
adding a particle or a hole to the Mott insulator phase, 
we show that 



G(T) = (V'(r,r)V'*(r,0)) 

-Fpr/fesT for r > 
-F,\r\/ksT forT<0 



(4.4) 
(4.5) 



in the r — > oo limit. This Green's function describes an 
extra vortex line created and destroyed at position r for a 
time r when r is positive. For r < 0, it represents a hole. 
We identify Fp with the (free) energy cost per unit length 
of a particle and Fh with the (free) energy cost per unit 
length of a hole in the Mott insulator. These energies 
will be lower than the energy per unit length of particles 
or holes localized to a single site {Ep and Eh), because 
they can gain entropy from wandering and hence lower 
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their free energies. In the quantum analogy, these 
energies are the "gaps" associated with particle and hole 
excitations in the Mott insulator phase. 

Th e Hubbard-Stratanovich transformed action of Eq. 
(4.1) immediately yields 



G(r) 



G(k,r), 



(4.6) 



with 

G(k,T) 

_ 1 



duj 



j-i(k) 



Eh+iujknT 



ra+1 



(4.7) 



E„-iukBT 



For r > 0, i.e., in the Mott insulator phase, the integrand 
has two poles, one each in the upper half and lower half 
of the complex-w plane. Thus both particles and holes 
are represented by the single complex field "ip- This is to 
be contrasted with the situation usually encountered in 
vortex physics, e.g. the one considered in Ref. where 
the coherent state complex field represents only one type 
of defect (vortex hues). 

Upon solving for the location of the poles and then car- 
rying out the momentum integrals by steepest descents 
in the limit |t| — )■ oo, we find that the defect free energies 
are 



(Afiif+EpEh 



Art 



sinh^(/i^/2) + sinh^{hy/2) 



cosh hj. + cosh h,. 



1/2 



TA/ii, (4.. 



where the upper sign applies to particles and the lower 
to holes, with 



and 



= fi- Ml, 



^ll ^ C/(n- 1/2) - 2t. 



(4.9) 
(4.10) 

(4.11) 



When /i > /zi, particles are preferred over holes, while 
for fi < ^i, holes are preferred. If we set fi = /^i, parti- 
cle/hole symmetry is restored. When h = Hj^ = 0, then 
as we approach the Mott insulator/supcrfluid transition 
(r — > 0^) along a generic path with /i 7^ /ii, the energy of 
the preferred defect vanishes like r, while the remaining 
defect energy remains finite as r ^ 0. When fi — fii, 
however, both defect energies vanish simultaneously pro- 
portional to ^/r. 

If we restore a transverse field, i.e. take h 7^ 0, then 
neither the particle nor the hole free energies defined by 
Eq. (l.J) vanish as r ^ 0. However, as discussed in 
Sec. [Ill C and illustrated in Fig. 5, we expect that tilted 



particles or holes proliferate at the transition at a nonzero 
angle relative to the z axis. To see that the corresponding 
defect energy vanishes at the transition, consider 



G(/,0,0) = (V(r,T)^*(O,O)) 
d'^kduj 



(2^)3 

exp[-F(0, 0)VfcBT] 



(4.12) 

,{k^x+kyy-u.r)Q^^^^^ (4.13) 

(4.14) 



where we parametrize the tr ansverse and r-axis separa- 
tions embodied in Eq. ( 4.12 ) as 



r 



{I sin 6 cos 0, 1 sin 9 sin t 
I cos 0, 



(4.15) 
(4.16) 



and we eventually take the I 00 limit. 

To carry out this Fourier transform, we rotate to a new 
coordinate frame (k'jCj'), with 



Lu COS 6 — krr sin 6 cos ( 



so that 



G(/,0,0) 



d^k'dw' 
(27r)3 ' 



ky sin Q sin < 



'G(k',w')- 



(4.17) 



(4.18) 



We note that the ui' integral will be dominated in the 
I — > 00 limit by the closest singularity to* of G(k, w) to 
the real axis in the lower half of the complex frequency 
plane, and will behave as exp[— F(k', 6*, 0)//fcBr], where 

F(k',9,(j)) :^iu*kBT. (4.19) 

The k' integrals are done by steepest descents, so 



(4.20) 



where k'^ is chosen such that dFCk' ,9,4>)/d'k'\^, = 0. 
F{9, (p) can be found numerically for any f?, 0, but the 
low energy excitations can be found quite simply. We 
look for solutions which obey 



F(6l,, 



as r 0^ 



(4.21) 



In other words, the lower-half plane pole of the Green's 
function should vanish as we approach the transition. 
But G~^(k' = 0,^^' = 0) = at the transition, and in- 
deed this is the only mode to vanish for r > 0. Thus we 
conclude that if (0, (j)) is chosen appropriately to provide 
a low energy excitation, then k'^ = 0. 

The angle that satisfies this condition is therefore de- 
fined by 



9F(k',6',( 



ak' 



= 0. 



(4.22) 



k'=0 



Hence, in a small wavelength expansion of G^^(k',w'), 
there must be no terms linear in k'. But notice that in 
the (k, uj) frame of Appendix ^ the action has no terms 
linear in k. Therefore, (k,a>) — (k',cj'), and the angles 
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(0, 4>) at which ther e are low e nergy excitations are ex- 
actly those in Eqs. ( |3.17 ) and ( 3.18| ). As expected, the 
direction in which there are low energy excitations is the 
same as the average direction of the defects in the super- 
fluid phase. 

To find the energy of these excitations, we need to find 
the poles of the Green's function at k' = and small 
Lu. From Eq. (CI), we immediately see that the pole of 
interest is at ilo' = r/[v\, so 



r 



(4.23) 



in the r — !■ 0"*" limit. 

This analysis is clearly valid away from the multicriti- 
cal point v = 0. If v = 0, then dF{\i! , 6, 0)/ak\,^o = 
for all 9 and 0, so there are low e nergy excitations at all 
angles. It is then clear from Eq. (|3.1l[) that, at the level 
of Gaussian fluctuations. 



for all and ( 
on the angle. 



F{9,(t>)^y/P (4.24) 
as r ^ 0^, with the coefficient depending 



2. Defect free energy at the XY transition 

While we expect the Gaussian fluctuation analysis to 
be valid (up to logarithmic corrections) near the generic 
transition, there will be nontrivial corrections near the 
XY transition at the multicritical point c^, = 0, cj^ = 0. 
Close to this special point, we find that the free energy 
FxY of defect lines at all angles scales in the same way 
with r within the Gaussian fluctuation analysis 



F 



XY 



(4.25) 



where ^xy ~ l/r is the Gaussian or mean fleld theory 
correlation length. To incorporate the effects of ther- 
mal fluctuations near the transition, we note that Fxv 
is actually an inverse correlation length in the imaginary 
time direction. Because the dynamic scaling exponent 
near the special multicritical XY transition is z = 1, we 
expect 



FxY ~ {r 



(4.26) 



vortex lines in the Mott insulator phase. Although the 
linear resistivity vanishes, supercurrents in the (x,y)- 
plane generate a nonzero voltage due to thermally ac- 
tivated nucleation and growth of particle-hole pairs. To 
calculate this resistivity, we estimate the free energy of 
a particle-hole pair of a given size in the presence of an 
externally imposed current J. We expect that on scales 
much larger than ag, the lowest energy configuration will 
be approximately given by a loop (see Fig. 7). Then, 
following the analysis of Fisher, Fisher, and Huse for a 
conventional Meissner phase, we estimate the free en- 
ergy of a loop of radius i? 3> oq lying in a plane normal 
to the current to be 



loop 



2ttRF 



(4.27) 



where the second term arises from work done against the 
Lorentz force and the first is the line free energy of the 
loop. If J is in the (0j = 7r/2,0j) direction, then the 
average energy per unit length of the particle-hole pair is 



F 



F 



1 

as illustrated in Fig. 7. 



Fie 



de, (4.28) 



(6=71/4, (^=0) 




(6=271/3, (^=0) 



particle 



(9=71/6, (t)=7l) 



FIG. 7. A particle/hole defect loop which gives rise to a 
nonlinear resistivity by nucleating and growing in response 
to a current perpendicular to the (a;, r) plane. 



where vxv « 2/3 is the correlation length exponent of 
the three-dimensional XY model, and is the transition 



point renormalized by fluctuations. Equation (4.26) is, 
of course, the behavior predicted for the energy gap near 
the transition predicted by Fisher et al. [R^ 



3. Nonlinear resistivity in the Mott insulator phase 

We now determine the form of the nonlinear current- 
voltage characteristics for a current perpendicular to the 



-Fj^Qpp has a maximum at ~ cF / J(t>o, yielding a free 

2 

energy barrier of approximately Trci^ /J4'o- Thus the 
voltage generated by loops going over the barrier takes 
the form 



with 



V 



Jt 



'Jt/J 



TTCF^ 



(4.29) 



(4.30) 
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We expect this to be valid for 
J < 



(4.31) 



Th erefore, near the XY transition, we find from Eq. 
([4.26D that 



Jt ^ [r - rc) 



(4.32) 



Near the generic transition, Jt does not vanish, but in- 
stead approaches a nonzero constant. The exact value 
of F is tedious to obtain, but it can be estimated as 
F^^iFp + Fh), or 



F^ I {Anif + EpEh 



Art 



^sinh^{h^/2) + smh^{hy/2) 
cosh + cosh h„ 



1/2 



(4.33) 



where A/ii is given in Eq. (4.9) 



B. Finite size effects in the Mott insulator 

In this section, we discuss the difference between the 
periodic boundary conditions appropriate to real quan- 
tum bosons and the free boundary conditions appropriate 
to flux lines. In the Mott insulator phase, the transverse 
magnetic field is screened in the bulk of the sample where 
the vortices are localized on lattice sites, a phenomenon 
known as the transverse Meissner effect. |l^ However, 
near the boundaries of the sample along the r-direction, 
the transverse field will be able to penetrate. We show 
here that surface effects fall off exponentially with dis- 
tance and calculate the corresponding screening length. 
This length scale plays a role similar to the usual London 
penetration depth in the Meissner phase, but is physi- 
cally distinct, as it arises from the existence of a periodic 
pinning potential. 

It is convenient to characterize surface effects by (7ii)r, 
where {'Hi)t is the matrix element of Eq. ( |3.3D in the 
Heisenberg representation. 



{Hi 



(4.34) 



and lipi) and \4>f) are initial and final boson states appro- 
priate to "free" boundary conditions. By "free" bound- 
ary conditions, we mean that one integrates freely over 
the positions at which vortices enter and exit the sample. 
As discussed in Ref. this boundary condition can be 
represented in the boson mapping by 



|V'.) = l^/>=eV^|0), 



(4.35) 



where no = J2j ^^'^ |0) is the vacuum state. 

Recall Hi is the hopping term, which leads to a free 
energy reduction due to the wandering of the flux lines. 



Because vortex fluctuations are less constrained near a 
free surface, we expect that Tii exceeds its bulk value 
near the surface r = 0. In Appendix we show more 
precisely that 



{'Hi)t - {Hl)r=o 



-r I T* 



as T — > CO, with 



Fh 



(4.36) 



(4.37) 



where Fp and F^ are the h-dependent gaps given by Eq. 
Equation (4.37) implies that the influence of a 



free surface extends over the distance a particle-hole pair 
(with free energy Fp + Fh) can exist in thermal equilib- 
rium at temperature T, as illustrated in Fig. 8. 




bote 



particle 



V 



FIG. 8. The screening length r* of a transverse magnetic 
field is determined by the distance over which a particle- 
hole pair can exist in thermal equilibrium. 

Away from the multicritical point, this particular 
"healing length" does not diverge near the superfluid 
phase: suppose, for example, that particles are the fa- 
vored defects whose energy vanishes as we approach the 
transition. However, any particle excitation near the 
boundary must be accompanied by a hole, with a much 
higher free energy, yielding a finite healing length. Once 
we cross the transition into the flux liquid or "superfluid" 
phase, the free energy per unit length of extra particles 
becomes negative. Only then can it extend across the 
whole length of the sample, and no longer requires an 
accompanying vacancy. 

For penetration at surfaces in other directions, we ex- 
pect (although we have not done explicit calculations) 
that the magnetic field falls off exponentially with a de- 
cay length of 



F{9,(j))+F{7T-0,(I> + Tr)' 



(4.38) 
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where {9, (p) is the angle of the normal to the surface. 
Again, this will not diverge unless Cr = = 0. 



V. SUPERFLUID PHASE 



To study the superfluid phase, we must include quar- 
tic powers of It/) in the action. Near the transition, the 
truncation (3.11) should be an adequate approximation, 
because other terms, such as and gradient operators 
combined with tp* ijj* ^ijj , are irrelevant variables. 

Note first that in the superfluid phase, both particles 
and holes are present in equilibrium, even away from 
the exceptional multicritical point. As illustrated in 
Fig. 9, {^{r,T)) ^ p-F+fknT^ where is the free en- 
ergy associated with creation of a vortex "head" or mag- 
netic monopole at position (r, r). Similarly, ('0*(r,r)) ~ 
q-F-Z^bT ^ where F_ is the free energy associated with 
creation of a vortex "tail" or magnetic antimonopole. |2^] 
In the Mott insulator phase, these heads and tails are ac- 
companied by long strings of (energetically unfavorable) 
particle or hole excitations, so and F_ diverge as the 
sample thickness L ^ oo. The presence of both particles 
and holes in the superfluid phase follows because and 
F- are finite, i.e.. 



\r\ 



lim (^(r,T)V'*(r,0)) = (V'(r, r))(iA*(r, 0)) ^ 0. (5.1) 




Hole 



<t|/(r,T)> 7^ 




Particle 
<t|/*(r,T)> 7i 



When c 7^ or 7^ 0, we find more precisely that 



FIG. 9. Vortex "head" and "tail" excitations (monopoles) 
associated with (-i/)) and (^*). The corresponding free en- 
ergies F+ and F- are finite in the superfluid phase, leading 
to (i/j) 7^ and (ip*) 0, while in the Mott insulator phase 
there is a free energy cost per unit lengtti, giving rise to 
{ip) = {ip*} = as Lz ~* oo. 



lim ('(/>(r,T)V'*(r,0)) = - |rlnr| (5.2) 



\r\ 



just above the transition. The explicit expression for 
the condensate fraction uq = is given in Ap- 

pendix ^ Thus both particles and holes will proliferate, 
even when only one of their energies vanishes at the tran- 
sition, similar to results found for vacancies and intersti- 
tials at the supersolid transition by Frey et al. for super- 
solids ||2l[], and as confirmed via the decoupling method 
of Appendix The ratio of particles to holes in the su- 
perfluid phase is more difficult to compute, and we rely 
on the decoupling method of Appendix ^ for a mean field 
estimate of this quantity. More generally, we expect that 
the probability of an m-fold particle or hole behaves as 



lim mx,,T)r[r{x^M 



(Irlnrl)", 



(5.3) 



which has logarithmic corrections to the mean field re- 
sults obtained in the decoupling approximation of Ap- 
pendix 



A. Magnetic field penetration 

Once in the superfluid phase, the magnetic field B 
deviates from its frozen Mott insulator value Bq — 
{n(/)o/aQ)z, unless the parameters are tuned to the mul- 
ticritical point. Let (5B = B — Bq. We show exphcitly 
that the excess magnetic field (5B is determined by the 
angle that which tilted particles and holes enter at the 
transition. 

The calculation of the magnetic field appears in Ap- 
pendix 1^. The excess field (to lowest order in the con- 
densate fraction) is given by 



(5B — (pov 



|r lnr| 
WttD 



(5.4) 



as r ^ 0, wher e v is defined in Eq. ( [3.16 ) and D is 
given in Eq. (C2). This adds further support to our claim 
that the transition is characterized by the penetration of 
"tilted defects," as the excess magnetic field is indeed 
in the direction of these putative defects. Note that the 
magnitude of this excess field has the same dependence 
on the distance from the Mott insulator phase as it does 
in the pure case near the Meissner phase at Hd ■ m 
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B. Excitation spectrum 



ACKNOWLEDGMENTS 



We now discuss the spectrum of excitations in the su- 
perfluid phase. In the superfluid or entangled flux Uquid 
phase, {tp{r, T)ip* (0, 0)) exhibits long range order, as was 
discussed at the beginning of this section. Here, we dis- 
cuss the flux lines themselves, and focus on correlations 
in the local vortex density 



(5.5) 



which are controlled by the excitation spectrum of the 
equivalent quantum superfluid. ^ Specifically, we expect 
that the spatial Fourier transform n(k, t) of the density 
has the correlation 

5(k,r) EE (n(k,r)n*(k,0)) (5.6) 
- 5(k,r = 0)exp[-eB(k)T/fcBT] (5.7) 

in the t —i- oo limit, where the decay in the imaginary 
time direction is controlled by the phonon-roton spec- 
trum eB(k) of the superfluid. [g^ Our interest here is to 
find the form of this spectrum when the magnetic field 
is tilted to make the problem non-Hermitian. We as- 
sume the tilt is large enough to destroy the transverse 
Meissner defect and apply the Bogoliubov approxima- 
tion 1^ to Eq. (O), which will be accurate deep in 
the superfluid phase. This approximation treats creation 
and annihilation operators in the zero momentum Fourier 
mode as c- numbers, and then neglects interactions be- 
tween bosons that are both outside the consdensate. Af- 
ter some straightforward calculation s, w e flnd that the 
excitation spectrum that enters Eq. (5.7) is 



eB(k) J2UnoeR{k) + e|(k) + iei{k), 



(5.8) 



where ng is the average vortex density, 
ejj(k) — 2<[cosh/i2,(l — cos/cj-oo) 

+ cosh /ij^(l — cos fcyflo)], (5.9) 

and 

e/(k) = 2t[smhhx sinkxaQ + sinh hy sin kyaa]. (5.10) 

The real part of the spectrum (^) reduces to the usual 
lattice Bogoliubov spectrum when h = 0. The imaginary 
part can be u nde rstood by considering the limit of small 
k, when Eq. (5 



can be written 
eB(k) = v/2C/noefl;(k) + iu^ ■ k, 



(5.11) 



where 



Uh = 2aot(sinh /i^E, sinh /ij,) (5-12) 

is parallel to th e tilt direction Ch defined in Eq. ( |3.13 ). 
Equation ( 5.12 ) is just the famous Landau formula for 
excitations in a moving superfluid, generalized to 
imaginary time. The imaginary part is required here to 
describe the drift of the entangled vortex lines relative 
to the column direction. This complex spectrum is sim- 
ilar to the spectrum of eigenvalues found for superfluid 
columnar in positions in the non- interacting case. |p^ , p7t . 
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APPENDIX A: MEAN FIELD DECOUPLING 



Following Sheshadri et al. Q we make the decoupling 
approximation 



alaj w {a]) a J + a]{aj) - (aj)(aj). 



(Al) 



so that the Hamiltonian ( p~l| ) takes the form 7t! = Ti^, 
with 



= ^n,{n, - 1) - ^ih, - t*[ra\ + - W], (A2) 



and where the order parameter is '0 = {o-i)- The effective 
hopping strength t* is given by 



t* — 2t{coshhx -t- coshhy 



(A3) 



Since the sites are now decoupled in Eq. (A2), we focus 
on a single site, and drop the index i. We expand the 
ground state wave function of the site |$) in the occupa- 
tion number basis: 



|$)=^/(n)|n). 



(A4) 



n=0 



where ^ |/("')P = 1 to insure proper normalization of 
1$), and |/(»T.)P is the probability of finding n particles 
on a site. Then 



^ = (a,) -^V^^r(n)/(n+l), 

oo 

r = (4) =^V^/(n)r(n+l), 



(A5) 
(A6) 



n=0 



and 



U 



E^{n)=Y.[-n{n-l)-^^n]\f{n)\ 



(A7) 
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For any fixed set of magnitudes we can maximize 

l^/^l, and hence minimize E, by equating all the phases of 
the f{n). So without loss of generality we assume that 
all f{n) are real in the ground state wave function. 

In the Mott insulator phase with n flux lines per site, 
we expect f{n') — 5n,n' , and it follows from Eq. ( A5) 
that the order parameter "0 = 0. In the superfluid phase, 
/(n') will in general be nonzero for all nonnegative inte- 
gers n'; however, near the transition, we expect defects 
with either n-l- 1 or n — 1 flux lines per sites to dominate. 
Hence, we neglect all other types of defects and write 



/(n') 



if n' = n + 1 
if n' = n-l 
vT^^TT^^ if n' = n 
otherwise, 



(A8) 



where a and (3 are variational parameters. If we reexpress 

a and /3 as 



a = e(l + uj), 
/3 = e(l-^), 



(A9) 
(AlO) 



then e is the average number of particles and holes (e <C 1 
near the transition) and uj is t he particle-hole asymme- 
try. Upon expanding Eq. (A7) to linear order in e, we 
obtain 



E^e[{Ep + Eh)+uj{Ep-Eh) 



(All) 



with Ep and Eh given by Eq. ( |3.14| ) and Eq. ( |3.15| ). By 
minimizing the bracketed term with respect to uj, and 
equating the result to zero, we find a transition from 
the Mott insulator to the superfluid in agreement with 
the Hubbard-Stratanovich result (LS). The particle- hole 
asymmetry is given by 



where 



/io = U{n- 1/2) -t*/2. 



(A12) 



(A13) 



We can also calculate the probability of finding more 
exotic defects, e.g., two extra particles on a site, as a 
check on the ansatz (A8). To this end, we minimize Eq. 
(A7) with respect to n' , with n' ^ n, subject to the con- 
straint 



(A14) 



and obtain 



_ t*^[f{n' + 1)V^7TT + f{n' - l)vn', ^^^^^ 



Cn' + t*%l)K 



for n' ^ with 



^ = ^ V^7Tlf{n')f{n' + 1), 



n'=0 



[/(n+l)V^rTT + /(n-l)V^] 



fin) 



(A16) 
(A17) 



and 



Cn' = ^("' - n)^ - - U{n - 1/2)] (n' - n). (A18) 

Although these equations are difficult to solve analyti- 
cally, it is possible to extract the power law dependence 
of f{n') on e near the transition. Because ~ k -y/e, 
we have a recursion relation 

f{n')^^e[f{n' + l) + f{n'-\)l (A19) 

with asymptotic solution 

f{n') - el"'-»l/2 (A20) 

as e — > 0. Thus the fraction of sites with n' particles per 
site falls off like ~ e'" ^"1 near the transition; the same 
result follows from the Hubbard-Stratanovich method. 



APPENDIX B: EXPANSION OF EFFECTIVE 
ACTION 



To expand the effective action of Eq. (3.6) in powers 
of tjj, first define the Fourier transform 



V'(k, w) = / dT^,{T)e 



-i(k-Xi —cut) 



and 



r{K^)^y2 r drrAr)e'^''"''-'^^^ 
. Jo 



(Bl) 



(B2) 



(B3) 



with periodic boundary conditions in all directions. 
Then, in the thermodynamic limit, the first term of Eq. 
.61) reads 



PJ2 f dT{J'%,rdr)^Ar) 



with 



J-i(k) = 



Pal / ^|0(k,c.)|^J-i(k), (B4) 



2t cos{kxao + ihx) -\- cos{kyao + ihy) ^ ^ 
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We shall also need the second term of Eq. expanded 
to quadratic order in ip^ namely 



X [V',(Ti)V*(r2)(r,aI(Ti)a,(r2))^ 
+ ,PUri)Mr2)(TradTi)al{T2) 



(B6) 



where we have used the fact that Tip conserves particle 
number. 

Upon changing variables so that the earlier time is r 
and the later time is r + St, Eq. (B6) becomes 



f-L/2 



/3^V / dT f ddT 

+ il;* (t + ST)^|;^{T) (a^{T + ST)al{T) 



(B7) 



The expectation values in Eq. (B7) can be eval uate d 
using the properties of the single site Hamiltonian (3.2), 
which is diagonal in the occupation number representa- 
tion. For example, 

a|(T + (5r)ai(T) 



E 

m— 



E' 

m=0 



(B8) 



where 



U 



m(m — 1) — /zm. 



(B9) 



In the limit L oo, the sums are dominated by m — n 
such that 



{n-^)<jj<n, 



and we have 

(^a,(r + (5T)a|(r)^^ = (n + l)e~ 

with Eh a nd i?p the hole and particle energies given by 
Eqs. (Ell) and (|]l|). 



(BIO) 

(Bll) 
(B12) 



We can now Four ier transform (B7), which upon com- 
bination with Eq . (B4) yields the action S{ip) that ap- 
pears in Eq. (3.5) to second order in ip: 



J-i(k)- 



(27r)3 
n 



n + 1 



Eh + iLoksT Ep — iuiksT 



(B13) 



This quadratic approximation to the action is used to 
study the Mott insulator phase in Sec. IV. To determine 
the universality class of the Mott insulator/superfluid 
transition, we need the long-wavelength contributions 
(up to quadratic order in k and to) to the coefficient of 
|?/'(k, w)p and the zero- frequency part of in the ex- 
pansion of Eq. (^^). It is tedious, but straightforward, 
to find the result used in Sec. HI, namely 



SW ~ Pal 



d'^kduj 



(27r)3 



|V'(k,aj)p[r - iCf^uj -f ich • k 



where 



1 



Dijhkj]+u j d^vdT\Tp{Y,T)\'^, (B14) 



1 



2t cosh hx -\- cosh 



n+1 
E^ 



n 



El 



n + 1 

En 



Ch 



flQ X sinh + y sinh hy 
2t [cosh/la; -I- cosh/ij,]2 



and 



(kBT)' 



n+1 

~E^ 

p 



El) 



The matrix Dij reads 



D 



Exx Exy 

Dxy Dyy 



with 



Dx - 



1 + cosh hx cosh hy ~ sinh 



(cosh hx + cosh hyY 



2t 



Aq 1 + cosh hx cosh hy — sinh^ hy 



2t 



(cosh hx + cosh hy)^ 



Oq —2 sinh hx sinh hy 
2t (cosh hx + cosh hy)^ ' 



(B15) 
(B16) 
(B17) 

(B18) 

(B19) 

(B20) 
(B21) 
(B22) 



while 



u = (3aQ 



(n + l)^ n^ n{n+l)iEp + Eh) 



E^ 
p 



El 



^p-^h 



{n + l)[n + i) n{n-l) 
EpE2p ElE2h 



(B23) 



where i?2p and E2h are the energies of "double particle" 
and "double hole" excitations respectively. 



E2h = 2^l - i2n - 3)U, 
E2p = {2n+l)U-2fi. 



(B24) 
(B25) 



Note that u > because E2h > 2Eh > and E2p > 
2Ep > 0. 
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APPENDIX C: RENORMALIZATION GROUP 

We first review tiie renormalization group treatment of 
the action of Eq. (3.11) away from the multicritical point 
c = Ch = 0. The analysis simpHfies if we transform the 
axes of space and time so that the imaginary time axis 
is no longer the z axis but coincides with the preferred 
direction of the particle or hole defects. We denote quan- 
tities in this rotated frame with a tilde. In this special 
coordinate frame, the term linear in k vanishes and the 
action takes the form 



S{ip)Kl3all \tP{\<i,Cj)\'^[r-iCCj + Dk^ 



k,Lj 

u I |V'(f,r)|^ 



(Cl) 



where C = |v| from Eq. ( [3.16 ). In general the quadratic 
wave vector dependencies will take the form Dijkikj. 
However, we can rotate the spatial axes so that the ma- 
trix Dij is diagonal, and then rescale the x and y axes 
(holding the area fixed) so that Dij = DSij, giving the 
simpler expression above, with 



2t\\r\{coshhx + coshhy)^ 

cosh ft,y)(cosh hx cosh hy 



X < — (cosh/i 



1) 



(cosh^ hx + cosh^ hy + 2) (cosh cosh hy ~ 2) 

1/2 

■ • (C2) 



(It can be shown that the quantity in braces is always 
positive when r > 0.) Also, there should appear an (D^ 
piece and an wk piece. However, these terms can easily be 
shown to be irrelevant variables near the transition, simi- 
lar to possible k and a)-dependent terms in the quadratic 
coupling. We add a high momentum cutoff, and thus 
allow only |k| < A. Following, e.g., Ref. we will in- 
tegrate out degrees of freedom ipCk^Lu) with momenta in 
a shell Ab~^ < |k| < A (and any uj), and rescale positions 
by b and imaginary times by 6^ to restore our original cut- 
off. We also rescale our fields by a factor 6'». In d = 2, u is 
marginally irrelevant at the Gaussian fixed point u — 0, 
allowing us to work perturbatively in u to lowest order. 
The most important diagrams for renormalizing r and u 
are shown in Fig. 10. In d -f 1 dimensions, we set z — 2 
and C — —d/2 to keep C and D fixed. Upon letting b — 
with I infinitesimal, we obtain the renormalization group 
recursion relations in = 2 



dr 
du 



(C3) 
(C4) 



Although u is irrelevant at the Gaussian fixed point, its 
slow decay to zero leads to logarithmic corrections to re- 
sults calculated in the mean field approximation. 




(a) 




(c) 

FIG. 10. (a) The most important diagram renormalizing 
r. (b) The most important diagram renormalizing u. (c) 
A diagram that vanishes identically. 

To exhibit these logarithmic corrections, we first solve 
the recursion relation for u to get [with uq = u{l = 0)] 



u{l) = 



Mo 



l + uol/{ATiCDI3'^aiy 
which behaves for large I like 



(C5) 



u{l) 



I 



Consider now the condensate fraction 
n, = mv,T))\^ 



(C6) 



(C7) 



right above the transition. We express no as a function 
of the parameters in the Hamiltonian, and iterate them 
under the renormalization group. Then 

7io(r, C, D,u,...)^ e^'^no{r{l), C{l),D{l), u{l), ...). 

(C8) 
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We iterate until r{l) « re^' — 1, i.e. I = l/2|lnr|, so 
that we are deep in the flux liquid or superfluid phase. 
Since uq « (3a'^r{l)/2u(l) far from the transition, we find 
using r{l) « re^\u{l) « AirCDP^al/l that 



\m 



no 



Irlnrl 



(C9) 



which has logarithmic corrections from the naive result 
no ~ r. 

We now turn to the case where = Ch = 0, so we 
are at the multicritic al po int. Up to a simple rescaling of 
lengths, the action (3.11) is then isotropic in space and 
time, the dynamical exponent is z = 1, and the univer- 
sality class is that of the three-dimensional XY model. 
As noted by Josephson , the superfluid density scales 
differently from the condensate fraction, with the con- 
densate fraction behaving as 



no 



{tc - r) 



2/3 



(CIO) 



near the transition while the superfluid density behaves 
as 



Ps ^ {tc - r) 



2/3-1/ »j 



(re 



\{d-2)i' 



(Cll) 



where rc — r is the distance from the transition. (The 
superfluid density is proportional to the inverse tilt mod- 
ulus in the flux line system. [ p^p5| ) We will have anoma- 
lous scaling for the superfluid density at the generic 
points of the phase boundary as well. Note that 77 = 
since we are at the upper critical dimension of the generic 
i?ci transition, so the corrections are limited to the log- 
arithmic terms. 



1. Validity of the continuum limit 

To explore the validity of taking the c ontinu um limit, 
we add a new term into the action of Eq. ( 3.11 ) that rep- 
resents a vestige of the periodic columnar array, namely 



where 



Jr(x) = p 



f2TTx\ 




cos 


+ cos 







(C12) 



(C13) 



We then investigate how the addition of this term changes 
the rcnormalization group. 

First we consider the case h = 0, so that we do not 
have to rotate our coordinate system as in Eq. (CI). 
Then, in momentum space, the additional terms all in- 
volve terms such as V'(k, a;)'0*(k -I- G,w), where G is a 
reciprocal lattice vector of the underlying square lattice. 
When we reach the point in our rcnormalization group 
proced ure w here (27r/ao)e' = A, all of the terms added 
in Eq. ( |C12 ) will be integrated out, leading only to finite 



changes in the remaining coupling constants. Hence, the 
critical phenomena at the Mott insulator to superfluid 
transition will be unchanged by the addition of the peri- 
odic potential. Note that this argument applies both to 
the multicritical point and to the generic transition with 
h = 0. 

When h 7^ 0, we must rotate our coordinate system as 



in Eq. (CI}. For most values of and Ch, all of the terms 



in Eq. (C12) will still involve a jump in the momentum 
k, so the argument above holds here as well. The excep- 
tion is when the direction of the defects is perpendicular 
to the columns, i.e., = 0. Then there will be terms in 
Eq. ( pl[ ) that involve a jump only in uj, not in k. Wc can 
rewrite these terms as 



55(V) 



/3ao2 



5-V'(k,u;) 



(27r)3 

X V*(k,(:' + A) V*(k,a) 



A) 



(C14) 



where the initial value of A is 2tt / a^. 

We now investigate how the parameters p and A be- 
have under the rcnormalization group. The diagrams 
renormalizing r, p, and u to one loop order are shown in 
Fig. 11. Note that all diagrams contributing to the rcnor- 
malization of p in this approximation vanish identically, 
as do all of the extra diagrams that renormalize r. This 
vanishing persists to all orders in perturbation theory, 
because any diagram that contains a loop (of more than 
one propagator) with all the arrows pointing in the same 
direction, as illustrated in Figs. 10(c), 11(a), and 11(b), 
must vanish identically. Moreover, all diagrams renor- 
malizing a propagator that contain any p-vertex must 
contain at least one such loop. Therefore, the zero-loop 
recursion relation for p is accurate to all orders in per- 
turbation theory, yielding: 



dp 

U 
dA 

la 



2p, 



= 2A. 



(C15) 
(C16) 



The recursion relation for r remains as given by Eq. (|C3|), 
accurate to one loop order. The recursion relation for u, 
however, has additional nonvanishing terms at one loop 
order, proportional to u^(p/A)^", for n a positive inte- 
ger. (There are also terms with extra powers of A in the 
denominator, but those vanish in the large I limit.) S ince 
p/A is independent of Z, only the constant in Eq. (C4) 
is changed; the structure of the equations remains the 
same. In other words, we expect that, aside from pref ac- 
tors, the results of our rcnormalization group calculations 
will be unchanged by the addition of a periodic potential. 
Therefore we can safely neglect terms with the periodic- 
ity of the columnar pins in our continuum description of 
the Mott insulator to superfluid transition. 
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If we now define a composite partition function 




(a) 




FIG. 11. (a) The diagrams involving factors of p that 
renormahze r. The dots are "interaction" vertices, with 



the interaction as given in Eq. (014). The diagrams all 
vanish identically, (b) The diagrams renormalizing p. 
These also vanish identically, (c) Nontrivial diagrams in- 
volving factors of p that renormalize u. 



APPENDIX D: FINITE SIZE EFFECT 

In this appendix we examine a particular type of free 
surface effect, by showing that 



{Hl)r - (Hi) 



(Dl) 



where {Tiijr is defined in Eq. ( 4.34 ). 

In the interaction representation, Eq. (4.34) takes the 
form 



(Hi>r = ^(T,exp 



-13 / dT"Hl{T') 



Mr) 



X Tt exp 



-/? / dr'Hi(T') 
Jo 



where 



(•>o = (V'/ 



-l3LHo 



,-PLHo 



and 







exp 


-/3 / dr'-Hiir') 




JQ 



(D2) 

(D3) 
(D4) 

(D5) 



X Tt exp 

then 



-/3 / dr'niir') 

J T2 

-13 [ ' dT'Hiir') 
Jo 



(D6) 



(D7) 



The Hubbard-Stratanovich transformation genera lize s 
straightforwardly to Z{ti,T2), and leads via Eq. (D7) 
to 



{H,)r = {K{t))s, 

where 

Kir)=Y,{J-%^*ir)^,ir) 



(D8) 



P / dr'Mr') 



Mr) 



5] I ^T, exp 

X Tr exp 

X (^r exp 
f^(^r') = Mr')al{r')+rAr')a.iT'), 



P I dr'Uir') 





P I dr'Mr') 





J / 



(D9) 



(DIO) 



>s - ^ y l[VMr)T^i^*ir) . exp[-5(^)], (DU) 



and 



SW^PY. r dT{J-%r^{r)^Ar) 



In ( Tt exp 



P / dT{Mr)al{T) 



+ r^{r)d^{T)} 



(D12) 



Alth ough Eq. ( D12 ) may appear to be identical to Eq. 
( |3.6[ ), the meaning of (•)o as well as the boundary con- 
ditions on ip have changed. For "boson" boundary con- 
ditions, ip was restricted to be periodic, i.e., i/jiir = 0) 
-04 (t — PL). Following Ref. we have inserted bound- 
ary conditions appropriate to vortices in a superconduct- 
ing slab. 
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In the Mott insulator phase, it is again sufficient to 
expand the action to quadratic order in ■0, yielding 



(D13) 



where 5*0 (V') is given by Eq. (4.1) and SSi^ip) is an extra 
piece due to the changed boundary conditions in a Mott 
insulator with n particles per site 



(2^)^ 



^(k,c.)^*(k,c.') 



^iuj' I3L 



[Eh + iLukBT){Eh 



iLu'ksT) 



[Ep — iu}kBT){Ep ~ iLo'kBT) 



(D14) 



We can now evaluate (7ii)r diagrammatically. Denote 
each SS{ip) with a cross. The diagrams that give the 
dominant contribution to (7ii)r — {'Hi)t=oo as r — > oo 
are those in Fig. 12, i.e., those with two external lines 
and an even (but nonzero) number of crosses. In the 
limit T — !■ oo we find 



{ni)r-{Hl)r-- 

which is Eq. ( |4.36| ). 



(D15) 



+ 



+ 



+ ... 

FIG. 12. The diagrams that give the dominant contribu- 
tion to {'Hi)t — ('Wi)oo as r ^ oo. The crosses are free 



boundary insertions as given in Eq. (D14) 



APPENDIX E: CALCULATION OF MAGNETIC 
FIELD IN SUPERFLUID PHASE 

To calculate the magnetic field in the superfluid phase, 
we use the thermodynamic relation 



B = -47r 



a/ 
9H' 



(El) 



where / is the free energy density / — 
— (ksT I N a^ L z) I nZ and the partition function Z is 
given in Eq. and Eq. (3.6). 

We use Eq. (3.5) to evaluate B. The piece that comes 
from Zq is easily seen to be exactly the magnetic field 



that we expect for a Mott insulator with occupation num- 
ber n. The overall field is 

4nkBT 



B 



— 



NalLz 



dS \ 



We first define 



and 



5B^B-% 



Then Eqs. (^ and (^ give 
6U 



00 V"o 



so that 



SB 



NalLz 



ao 



dS\ 



(E2) 



(E3) 



(E4) 



(E5) 



(E6) 



Near the transition, we can use the action of Eq. ( 3.11 ) 
to evaluate the averages in Eq. ( E6 ) , obtaining 



(SB 



^(i^(k,c.)r 



NLz J (27r)' 

4- 2ic'Luz — iDij{kiej + kjii) 



z J d2rdr^|V'(r,T)|*^ + •••, (E7) 



wher e the omitted terms reflect contributions left out of 
Eq. ( [3.1l| ). [Note that du/d^ is proportional to th e coe f- 
ficient of the uj\ip\'^ term that would appear in Eq. (3.11).] 
To see which terms will yield the most important contri- 
butions, we transform the axes of space and imaginary 
time as in Appendix We have already determined the 
scaling of the coefficients under the action of the renor- 
malization group, because they are the same coefficients 
as appear in the action. We discard those terms whose 
coefficients are irrelevant, and obtain 



6B 



"(2^ 



V'(k,cD) 



V - 2iD]i 



(E8) 



where v is given by Eq. ( 3.1q ). We then iterate the renor- 
malization group until \r{l)\ = 1, so that we are far from 
the transition, and can use mean field theory. If the 
flow begins on the Mott insulator side of the transition, 
then we will iterate until we are deep in the Mott in- 
sulator phase and Eq. (E8) vanishes. If the flow begins 
on the superfluid side of the transition, th en deep in the 
superfluid phase we can set |?A(f,f)| « ^y l3alr{l)/2u{l), 
independent of space and imaginary time. This yields 



6B « <^o/3aov(l^l' 



(E9) 
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We therefore identify v with the average direction of the 
tihed defects, and n = PaQ\v\{\ip\'^) with their density. 
Performing a scahng analysis similar to what was done 
for the condensate fraction tiq in Appendix we obtain 
(away from the multicritical point) 



SB 



|r Inrl 



Pov- 



(ElO) 



16nD 

using C — |v| and with D given in Eq. (|C2|). At the 



multicritical point, v = 0, so Eq. (E9) gives SB ~ near 
the transition. 
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